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Abstract. While numerically solving a problem initially formulated on an unbounded domain, one typically 
truncates this domain, which necessitates setting the artificial boundary conditions (ABC’s) at the newly formed 
external boundary. The issue of setting the ABC’s appears to be most significant in many areas of scientific comput- 
ing, for example, in problems originating from acoustics, electrodynamics, solid mechanics, and fluid dynamics. In 
particular, in computational fluid dynamics (where external problems present a wide class of practically important 
formulations) the proper treatment of external boundaries may have a profound impact on the overall quality and 
performance of numerical algorithms. 

Most of the currently used techniques for setting the ABC’s can basically be classified into two groups. The 
methods from the first group (global ABC’s) usually provide high accuracy and robustness of the numerical procedure 
but often appear to be fairly cumbersome and (computationally) expensive. The methods from the second group 
(local ABC’s) are, as a rule, algorithmically simple, numerically cheap, and geometrically universal; however, they 
usually lack accuracy of computations. In this paper we first present a survey and provide a comparative assessment 
of different existing methods for constructing the ABC’s. Then, we describe a relatively new ABC’s technique of ours 
and review the corresponding results. This new technique, in our opinion, is currently one of the most promising in 
the field. It enables one to construct such ABC’s that combine the advantages relevant to the two aforementioned 
classes of existing methods. 

Our approach is based on application of the difference potentials method attributable to V. S. Ryaben’kii. This 
approach allows us to obtain highly accurate ABC’s in the form of certain (nonlocal) boundary operator equations. 
The operators involved are analogous to the pseudodifferential boundary projections first introduced by A. P. Calderon 
and then also studied by R. T. Seeley. The apparatus of the boundary pseudodifferential equations, which has formerly 
been used mostly in the qualitative theory of integral equations and PDE’s, is now effectively employed for developing 
numerical methods in the different fields of scientific computing. 

Key words, infinite-domain problems, artificial boundary conditions, pseudodifferential equations, difference 
potentials method, auxiliary problem, boundary equations with projections. 
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1. Introduction. Artificial boundary conditions (ABC’s) furnish a widely used approach for 
the numerical treatment of boundary- value problems initially formulated on unbounded domains. 
These boundary conditions are typically set at the external boundary of a finite computational 
domain once the latter is obtained from the original unbounded domain by means of truncation. 
Implementation of the ABC’s enables completion of the “truncated problem” and, therefore, makes 
this problem available for solution on the computer. 

For almost any problem formulated on an unbounded domain, there are, generally speaking, 
many different ways of closing its truncated counterpart. In other words, the choice of the ABC’s 
is never unique. Clearly, the minimal necessary requirement of ABC’s is to ensure the solvability 
of truncated problem. If, however, we restrict ourselves to this requirement only, then we cannot 
guarantee that the solution found inside the computational domain will be anywhere close to the 
corresponding fragment of the solution to the original (infinite-domain) problem. Therefore, we 
must additionally require of the ABC’s that the two solutions be in a certain sense close to each 
other on the truncated domain. An ideal case here would obviously be an exact coincidence of 
these two solutions, which leads us to formulating the concept of exact ABC's. Namely, we will 
refer to the ABC’s as being exact if one can complement the solution calculated inside the finite 
computational domain to its infinite exterior so that the original problem is solved. The concept 
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of exact ABC’s appears to be useful for the theoretical analysis of infinite-domain problems. 

To provide a simple one- dimensional (ID) example of exact ABC’s, we consider a half-line 
problem with the compactly supported right-hand side (RHS): 


(1.1a) 


(1.1b) 

(1.1c) 


d 2 u o 

= *>0, 
supp f(x) C [0, X), 
u(0) = 0, 

u{x) — * 0, as x — ► +oo. 


Equation (1.1a) is homogeneous for x > X; therefore, it has (for x > X) two eigensolutions. 
The first eigensolution vanishes as x — • +oo and the second one infinitely grows as x — * 4-00. 
Boundary condition (1.1c) can be met iff the increasing mode (eigensolution) u^^x) = e^ x does 
not contribute to the solution of (1.1a) on the entire semi-infinite interval [X, +00). To prohibit 
this growing mode and to admit only the decaying one, u( 2 )(x) = e~^^ x , we use the following 
first-order homogeneous differential relation: 


( 1 . 2 ) 


du 

dx 


x=X 


+ \^\u 


:=X 


= 0 , 


which obviously provides us a desirable exact ABC at the artificial boundary x = X . Let us 
emphasize that relation (1.2) exactly transfers boundary condition (1.1c) from infinity to the finite 
boundary x = X. In other words, (1.2) adequately takes into account the structure of the solution 
to (1.1) from outside the finite interval [0, X] without explicitly calculating this solution. 

Although the above example provides some understanding of exact ABC’s, it is still not com- 
prehensive because of its ID nature. Therefore, let us consider the following two-dimensional (2D) 
problem for the Poisson equation (see also [1]): 


(1.3a) 


(1.3b) 


1 d ( du\ 1 d 2 u 

rlh\lb) + r*d(P~^ r,e ^ T ~ 0, 0 ^ e<2n ’ 

supp /(r, 0) C B C {r < i? 0 }, / f(r,0)rdrd0 = 0, 

J i B 

u(r, 6 ) — ► 0, as r — >+00. 


Again, equation (1.3a) is homogeneous for r > R 0 . Fourier transforming (1.3a) with respect to 0 
for r > Rq, we obtain 


1 d f duir\ k 2 „ „ 

(1 ' 4) ( r rf7)-^“* = 0 ' riR °’ * = o.±i.±2,... 

To satisfy (1.3b), we must prohibit the increasing eigensolution of (1.4), uj^(r) = t-W, on the 

entire semi-infinite interval [i? 0 , +00) for k = ±1,±2, . . . and allow there only the decaying mode 

u^\r) = For k = 0, we must prohibit both eigensolutions of (1.4), u l 1 * = const and 

* ( 2 ) 

Uq = lnr. Therefore, the following countable set of relations 

IJfcl . 

+ — Uk =0, fc = 0,±l,±2,..., 

T — Rq r t — Rq 


( 1 . 5 a) 


dv,k 

dr 
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(1.5b) ii 0 (R 0 ) = 0, 

constitutes here the exact ABC’s at the artificial boundary r = R 0 . One can easily make sure that 
boundary conditions (1.5) are spatially nonlocal (in the original variables). Indeed, since (1.5a) 
contains an absolute value of the wavenumber k : an inverse Fourier transform of (1.5) yields a 
nonlocal expression with some pseudodifferential operator ($DO). 

We emphasize here that the situation illustrated in the last example is fairly general. For many 
different problems, including those that originate from applications (e.g., from continuous media 
mechanics or from electrodynamics), the exact ABC’s are nonlocal (except in some trivial, e.g., ID, 
cases). In fact, the nonlocal character is one of the most essential features of exact ABC’s. Another 
essential feature is that in a vast majority of cases such boundary conditions can be derived easily 
only for some particular geometries. Indeed, returning to the second example above we see that if 
the shape of artificial boundary were not circular, then the direct implementation of the Fourier 
transform would be impossible. 

We will now briefly survey some work on ABC’s conducted by different authors over the recent 
years. First, we will concentrate on efforts to construct and numerically implement the exact 
ABC’s. As one will easily see, the main tool for constructing such boundary conditions is integral 
transformations (e.g., Fourier and/or Laplace); in so doing, the artificial boundary should have 
some regular shape (e.g., linear, circular, elliptic, etc.). It is also important to mention that this 
group of methods applies, generally speaking, only to linear problems. Therefore, while discussing 
the exact boundary conditions we mean exactness for the linear formulation. If, however, the 
original problem is nonlinear, then one can often linearize it under certain conditions. For example, 
linearization in the far field is a common approach in fluid dynamics. If this is really the case, we 
will retain the term “exact ABC”, assuming exactness within the accuracy of linearization. 

In [2, 3, 4], Engquist and Majda develop the ABC’s for some time-dependent problems, in 
particular, for wave propagation problems (wave equation, first-order hyperbolic systems) and 
for linearized potential transonic flows. Their approach is based on representing the solution as a 
superposition of waves and excluding all incoming waves from the solution at the artificial boundary. 
This is done exactly, using an apparatus of integral (Fourier) transformations. Note, ABC’s that 
employ an idea of prohibiting the incoming waves near the artificial boundary (such waves may also 
be interpreted as reflected ones) are frequently called non-reflecting boundary conditions (NRBC’s). 
The ABC’s [2, 3, 4] are based on analyzing the dispersion relations for the waves that contribute to 
the solution and on constructing special dispersion relations that correspond only to the outgoing 
waves to be solely left in the solution (as in the second example above, we select only the decaying 
modes in transformed space). Since these special one-way dispersion relations contain non-integer 
powers, ABC’s [2, 3, 4] in the physical variables are formulated using \PDO’s, i.e., they appear to 
be nonlocal in both space and time. 

Gustafsson [5] analyzes another hyperbolic problem that is more complicated from the stand- 
point of constructing the ABC’s. Namely, the support of initial data is no longer concentrated 
inside the computational domain (as in [2, 3, 4]) but is also permitted to spread beyond the ar- 
tificial boundary. Implementation of the Laplace transform in time and the Fourier transform in 
space yields in this case a nonhomogeneous (unlike [2, 3]) nonlocal exact ABC at the plane artificial 
boundary. Note, both boundary conditions in [2, 3, 4] and [5] are then approximated by some local 
relations for more computational convenience (see below). 

Sofronov [6] considers an example of noticeable interest for the wave equation; namely, he 
constructs the exact ABC’s in three dimensions at the spherical artificial boundary. The approach 
of [6] is based on using the Laplace transform in time and expanding the solution in space in 
terms of spherical functions. For the finite-difference formulation, standard spherical functions 
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are substituted by a special finite- dimensional orthonormal basis called the difference spherical 
functions [7]. 

Note, all above-mentioned techniques for the wave equation (see [2, 3, 4, 5, 6]) have strict 
limitations on the shape of the artificial boundary; although, the author of [6], for example, proposes 
some recipes for possible weakening of these requirements. 

Hagstrom and Keller [8] suggest to consider the exact ABC’s as a characterization of data 
at the artificial boundary in terms of belonging to certain admissible subspaces. The latter are 
specially defined to ensure the solvability of the exterior problem (i.e., the one formulated outside 
the computational domain) in an initially prescribed class of functions, e.g., bounded or vanishing at 
infinity. (Note, both the ID and 2D examples above provide this type of classification by selecting 
only the decaying modes as x — * oo or r — * oo; see (1.2) and (1.5), respectively.) For some 
problems formulated on cylinders, Hagstrom and Keller calculate these admissible subspaces on 
the plane artificial boundary normal to the cylinder element. They assume that the coefficients 
of partial differential equations (PDE’s) in the far field depend only on the transversal but not 
the longitudinal coordinate and use the separation of variables based on expanding the solution in 
terms of transversal eigenfunctions. The latter appear to be the Fourier harmonics for the simplest 
case of constant coefficients. In later work, Hagstrom (9, 10] and Hagstrom and Keller [11] further 
develop the technique introduced in [8] and extend its applicability. 

Givoli and Keller [12] propose the construction of (nonlocal) exact ABC’s for the Laplace 
equation and for some problems in elasticity. Their boundary conditions are based on the so-called 
Dirichlet-to-Neumann maps, which express the normal derivative of the solution at the artificial 
boundary in terms of boundary values of the solution itself. These maps (which sometimes are also 
called the Poincare- St eklov operators, see, e.g., [13]) are calculated in [12] analytically (using Fourier 
transform) for the circular and spherical artificial boundaries. In [14], Givoli and Vigdergauz use an 
apparatus of the Dirichlet-to-Neumann maps to construct the ABC’s for the Helmholtz equation 
and for the elastostatics system; the artificial boundary in [14] is composed of a semi-circle and 
two semi-infinite straight lines, which is a typical geometric setup for geophysics. In [15], an 
approach based on calculation of the Dirichlet-to-Neumann maps was extended by Givoli to treat 
time-dependent problems. The method of [15] is based on the analysis of the “stationary” system 
that arises on the upper time level when one integrates the time-dependent problem by a certain 
implicit method. Boundary conditions [15] are nonlocal (exact) in space and local in time. Finally, 
Givoli and Cohen [16] analyze an essentially time- dependent problem for the wave equation and 
for the elastodynamics system. From the standpoint of constructing exact ABC’s, time-dependent 
problems are generally much more difficult to handle than stationary ones. Indeed, in the time- 
dependent case the exact ABC’s generally appear to be nonlocal in both space and time (see, e.g., 
[2, 3]). This nonlocality may cause severe computational problems, mainly due to the constantly 
increasing amount of past information the numerical method must store in memory as the solution 
is advanced in time. However, for particular classes of problems, e.g., for hyperbolic systems with 
an odd number of space dimensions, one can use the existence of the lacunas (some explanation of 
this concept is given in Section 3) and limit the required amount of past information by some fixed 
value. This circumstance was effectively used by the authors of [16]. Moreover, since the Kirchoff 
integrals are used in [16] for constructing the ABC’s, those geometric restrictions (on the shape of 
artificial boundary) that usually apply in practical calculation of the nonlocal boundary conditions 
(see, e.g., [2, 3, 4, 5, 6]) are not encountered. On the other hand, calculation of the Kirchoff integrals 
has its own limitations, namely, it requires an explicit knowledge of the fundamental solutions. We 
also note that the possibility of using the lacunas to effectively reduce the required computational 
effort when calculating the exact ABC’s for time-dependent problems was earlier pointed out by 
Ryaben’kii [17]. 

Gustafsson [18] and Ferm and Gustafsson [19] investigate an inviscid flow of a compressible gas 
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in a duct with parallel walls. Linearizing the stationary Euler equations in the far field and Fourier 
transforming these equations in the cross-stream direction, they obtain an exact nonlocal ABC 
at the straight outflow boundary normal to the duct. The ABC’s [19] use the integral principle 
of conservation of mass and ensure the downstream boundedness of the solution up to infinity. 
The authors of [19] also justify the well-posedness of the corresponding problem if the boundary 
conditions [19] are directly used to treat the time-dependent case. In [20], Ferm modifies the above 
technique and obtains nonlocal ABC’s for both inflow and outflow boundaries in the channel; these 
boundaries are again the segments of a straight line. Moreover, an Engquist-Halpern approach [21] 
is implemented in [20] to accelerate convergence of the pseudo-time iterations to a steady state. 
The technique of [21] by Engquist and Halpern generally provides far-field boundary conditions for 
hyperbolic PDE’s. The idea is to combine exact ABC’s for the spatial part of the corresponding 
differential operator with some NRBC’s that make the artificial boundary transparent for outgoing 
waves. Numerical experiments by Ferm [20] corroborate that the Engquist-Halpern technique of 

[21] implemented on the basis of exact ABC’s for the linearized stationary Euler equations can 
significantly improve the convergence rate of pseudo- time iterations for duct flow problems. In [22], 
Ferm proposes an analogous (to [20]) approach to construct the exact ABC’s at the elliptic artificial 
boundary for inviscid compressible external flow problems. In [23], he modifies the approach of 

[22] to accelerate the convergence of pseudo-time iterations to steady state; since it turns out that 
the Engquist-Halpern technique of [21] is not as effective for external flows [23] as it is for flows in 
the duct [20], the modification proposed in [23] differs from that from [20] and is based on a slight 
perturbation of the free-stream Mach number. 

As mentioned in [20, 23], when solving a steady-state problem by pseudo-time iterations the 
direct im plementation of stationary exact ABC’s may result in a relatively slow convergence. There- 
fore, special acceleration techniques (see [20, 23]) are required to obtain an algorithm that would at 
the same time be highly accurate (exact ABC’s) and computationally effective (rapid convergence). 
However, it is very interesting to mention that when implemented along with a multigrid iteration 
procedure for the Euler equations (see Ferm [24]), nonlocal ABC’s [22] no longer slow down the con- 
vergence and therefore, do not require any additional acceleration technique (as in [23]). According 
to Ferm [24], the number of multigrid cycles required for reducing the initial error by a prescribed 
factor appears to be approximately the same for both ABC’s [22] and some local characteristic 
NRBC’s (to be described below). Our own computational experience [1] supports and expands 
on the above-mentioned results. Namely, a direct combined implementation of stationary exact 
nonlocal ABC’s [25] (see also [1]) and the multigrid iteration procedure [26, 27, 28] for integrating 
the Navier-Stokes equations results in a drastic convergence acceleration (sometimes by a factor of 
3) in comparison with some standard local NRBC’s. Boundary conditions [25] are constructed on 
the basis of the difference potentials method (DPM) proposed by Ryaben’kii [29, 30]; sections 2 
and 3 below are devoted entirely to the analysis of this approach. 

In [31], Verhoff, Stookesberry, and Agrawal construct ABC’s for inviscid compressible external 
flow computations. An interesting feature of this approach is that the Euler equations are linearized 
in the far field against the constant-pressure background; however, a special change of variables 
allows the nonlinear thermodynamic relations to be retained. This enables one to explicitly take 
into account entropy-wake solutions (i.e., rotational effects) that are relevant to inviscid treatment 
of the far field. Again, the Fourier transform (combined with a certain iteration technique) is used 
to solve the far-field equations and to obtain the ABC’s at the C-type artificial boundary that is 
composed of parabolic (inflow) and linear (outflow) segments. In [32], Verhoff and Stookesberry 
extend the above approach to duct problems, and in [33] Verhoff uses an analogous technique to 
treat O-type configurations for circular artificial boundaries. 

Among other papers devoted to constructing the (nonlocal) exact ABC’s, we mention work by 
Fix and Marin [34], in which the authors solve the Helmholtz equation in an axially symmetric 
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duct and construct the exact ABC’s on its lateral boundary expanding the solution in terms of 
the trigonometric and Hankel functions; work by Zorumski, Watson, and Hodge [35], and work by 
Hodge, Zorumski, and Watson [36], in which the exact ABC’s at the plane transversal artificial 
boundary are developed (again, using the separation of variables) for the Helmholtz equation in a 
semi-infinite rectangular constant-section duct; work by Watson and Zorumski [37], in which the 
method of [35, 36] is extended to treat ID time-periodic duct acoustic phenomena described by the 
linearized Euler equations in the far field; and work by Jiang and Wong [38], in which an exact 
NRBC containing $D0 is obtained at the plane artificial boundary for a general second-order 
hyperbolic equation, provided that the corresponding dispersion relation is known (the technique 
in [38] is analogous to that in [2, 3]). Nataf [39] calculates an external incompressible viscous flow 
past a finite body using the vorticity-stream function formulation of the Navier-Stokes equations. 
The linearized equations are solved outside the rectangular computational domain with the help of 
separation of variables. The resulting solution is then used to construct the exact nonlocal ABC’s 
for the Poisson equation that describes the far-field stream function. The boundary conditions for 
vorticity used in [39] are local; they are proposed by Halpern in [40]. Note, the far-field vorticity is 
described by the linear ad vection- diffusion equation; the exact ABC’s for this equation are derived 
in [40] using the cross-stream Fourier transform; then these boundary conditions are approximated 
by some local relations (see below). 

Clearly, the main advantage of using the exact ABC’s for computations is the high accuracy 
of the results. This statement is corroborated in many of the previously mentioned papers and 
generally remains true even when the exactness is regarded only within the accuracy of linearization. 
However, in reality the exact ABC’s can be achieved in numerical practice only rarely. A few 
reasons account for that. First, as could be seen from most of the papers cited above, the main 
apparatus used for obtaining the exact ABC’s is the integral transforms, which severely limits the 
class of admissible artificial boundaries. In practice, such geometric limitations often appear to be 
unacceptable since the artificial boundaries are typically defined by the discretization, i.e., by the 
grid, used for each specific case. As a rule, the grid is generated to reflect some essential geometric 
elements of the specific problem, for example, it can be fitted to the (inner) solid boundary. In 
so doing, the shape of the (external) artificial boundary may be rather complicated and cannot 
be easily modified for the convenience of setting the ABC’s. Second, even for those cases in 
which the exact ABC’s can be theoretically obtained, their practical implementation may encounter 
serious difficulties, mostly because of the considerable amount of computer resources (CPU time and 
memory) associated with the calculation of integral transforms. The corresponding requirements 
may still be reasonable for the steady-state problems but usually become too high for the time- 
dependent ones. 

These arguments justify numerous attempts by different authors to construct approximate 
local ABC’s. One widely used approach is simply to develop local approximations to the previously 
derived exact ABC’s. As mentioned above, the nonlocality is usually caused by the fact that exact 
boundary conditions in the transformed space contain some special expressions, such as absolute 
values (see (1.5a)) or radicals (the latter are relevant to the one-way dispersion relations in the 
wave propagation problems). Basically, each one of these special expressions is a symbol of the 
$DO that arises after the inverse transform. If one develops some rational (e.g., Taylor or Pade) 
approximation to such a symbol, then the resulting boundary condition in the physical variables 
acquires the form of a certain differential relation, i.e., it becomes local. The idea of constructing 
rational approximations to the symbols of ’J'DO’s was implemented by many authors, e.g., Engquist 
and Majda [2, 3, 4], Gustafsson [5], Jiang and Wong [38], Halpern [40]; as well as Jin and Braza 
[41] for the incompressible shear flows calculation, Blaschak and Kriegsmann [42] for the Maxwell 
equations (more precisely, for the Maxwell equations inside the computational domain coupled with 
the wave equation in its exterior), Kroner [43] for the linearized Euler equations, and Johnsen and 
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Lynch [44] for the shallow water problems. Note, in the work of Johnsen and Lynch [44], the ABC’s 
are obtained on the basis of the Klein-Gordon equation, which describes the dispersive waves; the 
case of dispersive waves is also analyzed by Higdon in [45]. 

Obviously, the ABC’s obtained using rational approximations to the symbols of $DO’s have 
the same geometric restrictions as the original nonlocal ABC’s. However, the new feature, namely, 
the locality of approximate ABC’s, presents a major advantage for numerics. 

The difference between nonlocal ABC’s and their local approximations can be considered in 
several different ways. Of course, the first question arising here is the question of convergence. 
A rigorous analysis of the convergence of local approximations to the corresponding fDO’s is 
provided by Hagstrom [46]. In particular, he shows that the local approximations used by Engquist 
and Majda [2] do converge to the corresponding ^DO’s over finite time intervals. Another useful 
approach is proposed by Higdon [47] for wave propagation problems. Namely, he shows that certain 
rational approximations to the corresponding $DO’s can be factorized so that the resulting local 
ABC appears to be exact only for the waves with some specific incident angles (rather than for 
all waves, as in the case of exact ABC’s). Higdon’s factorization result was further generalized 
by Jiang and Wong [38]. We should also mention that the ABC’s by Higdon [47] are constructed 
directly for the discrete formulation of the problem, on the basis of the dispersion relation for the 
finite- difference scheme. This eliminates one step from the general numerical procedure, namely, 
the need to discretize the continuous boundary conditions. 

An alternative approach to approximating the exact ABC’s consists of retaining only a few 
leading terms in the far-held asymptotic expansion of the solution and then using the obtained 
truncated expansion to set the ABC’s. This technique may essentially reduce the required compu- 
tational effort in comparison with the cost of the original exact ABC’s. The idea of the above type 
was employed by Sa and Chang [48] to set the ABC’s for vorticity when integrating the incompress- 
ible Navier-Stokes equations around a cylinder. Burkhart [49] and Burkhart et al. [50] derive an 
asymptotic expansion for the finite-difference fundamental solution of the three-dimensional Laplace 
operator on a Cartesian grid and then use a few leading terms of this expansion to set the ABC’s 
for an external flow problem that is solved within the full-potential framework. Wubs, Boerstoel, 
and Van der Wees [51] use a Fourier representation of the far-field solution to the Laplace equation 
(see (1.4)) to calculate a potential flow around an airfoil. The ABC’s [51] are again derived from 
the first few leading terms of the expansion; as the artificial boundary approaches the airfoil, more 
terms are required to maintain the accuracy. We also mention here the earlier work of Thomas and 
Salas [52], in which only one leading term of the far-field potential representation is retained; this 
term corresponds to the point-vortex model. 

Except for approximating the exact ABC’s, there are, of course, independent techniques for 
constructing the approximate local ABC’s. In particular, the approach introduced by Bayliss and 
Turkel [53, 54, 55] and Bayliss, Gunzburger, and Turkel [56] is also based on using the far-field 
asymptotic expansion of the solution. However, the authors do not directly use this expansion to 
set the ABC’s but construct special local differential relations that annihilate a certain number 
of leading terms in the aforementioned expansion. Being applied at the artificial boundary, these 
relations provide some local ABC’s. It is interesting to mention that sometimes local ABC’s [53, 54, 
55, 56] may coincide with those obtained by means of rational approximation to the ^DO symbols. 
We also note that an apparatus of asymptotic expansions for constructing the approximate ABC’s 
was extensively used by Hagstrom [9, 10], Hagstrom and Keller [11], Hagstrom and Hariharan [57], 
and Hagstrom [58]. 

Basically, both above-mentioned approaches to constructing local ABC’s (the one based on 
rational approximation to the 4>DO symbols and the one based on asymptotic expansion of the 
solution) provide an essential simplification of the numerical algorithm compared with the direct 
implementation of nonlocal ABC’s, as well as a substantial reduction in the required computer effort. 
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As concerns the accuracy of computations, it generally improves as the approximate conditions 
approach the exact ones. However, one restriction, namely, the requirement that the artificial 
boundary be of some simple shape, still holds here. Moreover, many authors mention that usage 
of the approximate local ABC’s may poorly affect well-posedness of the problem. In particular, 
certain types of rational approximations to the $D0 symbols may lead to instabilities, whereas 
some other approximations provide a stable truncated problem. Engquist and Majda [3] give some 
examples of both stable and unstable local ABC’s. Trefethen and Halpern [59] investigate and 
compare different ways of approximating the symbol of the ^DO by a rational function (Pade, 
Chebyshev, least-squares) and characterize stability in terms of algebraic properties of the rational 
approximant (locations and multiplicities of poles and zeros). For the case of hyperbolic systems 
with non-zero initial data outside the computational domain, the question of well-posedness of the 
truncated problem was studied by Gustafsson in [5] and more extensively in [60], 

Before proceeding to the description of some purely local approximate approaches, let us men- 
tion here the so-called perfectly matched layer technique [61] by Berenger. In a sense, this technique 
occupies an intermediate position between local and nonlocal methods. The idea is to surround the 
computational domain by a finite-thickness layer of a special model medium with such properties 
that the outgoing electromagnetic wave rapidly attenuates in this layer. The aforementioned spe- 
cial medium is matched to a vacuum so that the boundary produces no reflection for any incident 
wave. On one hand, the approach of [61] is local since it does not include any global relations along 
the boundary; on the other hand, it is nonlocal since it requires enlargement of the computational 
domain. (In this sense, the treatment of the far field proposed in [61] cannot be exactly called the 
boundary conditions.) Geometrically, the shape of the computational domain studied by Berenger 
is rectangular, with the sides of the rectangle aligned in the Cartesian directions. In practical terms, 
the technique of [61] demonstrates its superiority for electromagnetic waves calculation over the 
local boundary conditions [42] based on the Pade approximation. 

We now describe another approach for constructing the approximate local ABC’s. Because of 
its algorithmic simplicity, low computational cost, and geometric universality, this approach has 
become very popular in computational practice. At present, it is most widely used in different 
fields, in particular, computational fluid dynamics (CFD). 

Clearly, numerous nonstationary problems in CFD can be successfully modeled using the wave 
propagation “language”. As mentioned above, many authors that investigate unsteady problems 
of the wave propagation type suggest the construction of ABC’s that are based on the principle of 
allowing any wave that travels toward the boundary to leave the computational domain without a 
reflection. Being implemented with the help of integral transformations, this principle yields the 
exact nonlocal ABC’s (for those problems, in which the exterior solution is composed of outgoing 
waves only, see above). For time-dependent flow calculations, the idea of treating the external 
boundary remains the same as for general wave propagation problems: to make the boundary 
transparent for outgoing waves and to eliminate the spurious reflections of waves from the boundary 
back into the computational domain. However, a completely different way of implementation of 
this principle, which is based on the (quasi-)one-dimensional characteristic analysis, leads here to 
essentially local ABC’s. 

Hedstrom [62] considers the one-dimensional Euler equations (gasdynamic system) on a finite 
domain (segment). Using a Riemann representation near the boundary points, he explicitly calcu- 
lates the number of characteristics entering the domain and the number of characteristics leaving 
the domain. Clearly, characteristics that enter the computational domain correspond to incoming 
waves since the data are transferred inward along these characteristics from outside the domain. 
Analogously, characteristics that leave the computational domain correspond to outgoing waves 
since the data are transferred outward along these characteristics from inside the domain. In terms 
of setting the boundary conditions, the values of outgoing Riemann invariants should be extrap- 
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dated along the characteristics from inside the domain to its boundary. It makes the boundary 
transparent for outgoing waves, and in the special case of supersonic outflow (when no incoming 
characteristics exist), it enables one to specify all necessary boundary conditions. The situation is 
different for the subsonic case, when for both inflow and outflow types of the boundary there are 
always incoming, as well as outgoing, characteristics. Obviously, specification of only the outgo- 
ing data would leave the problem subdefinite; therefore, the incoming Riemann variables must be 
prescribed as well. This step, generally speaking, can be accomplished in a variety of ways and 
here one can most clearly see the difference between “simply the appropriate boundary conditions” 
and “the true ABC’s”. Namely, to simply obtain a well-posed problem on a finite interval, one 
can select the boundary conditions, i.e., specify the incoming Riemann invariants, from a fairly 
wide class of data (this problem has been studied by many authors). On the other hand, to con- 
struct the true ABC’s one has to take into account the structure of the solution from outside the 
computational domain. In particular, the exterior solution should be used to prescribe the values 
of incoming Riemann variables. The simplest and most commonly used approach here is to set 
all the incoming Riemann variables outside the computational domain constant and equal to the 
corresponding free-stream values. That is exactly what was done by Hedstrom in [62]. We should 
mention that this approach is fully relevant to one- dimensional problems only. The corresponding 
boundary conditions are sometimes called radiation boundary conditions. Hedstrom [62] shows 
that these conditions provide zero reflection for the rarefaction waves, while for weak shocks, the 
reflection is cubically small. 

The local radiation boundary conditions can also be applied to multidimensional problems; 
although some simplifying assumptions have to be done in regard to the behavior of the solution 
outside the computational domain. For example, for the case of two-dimensional Euler equations 
in subsonic regime one cannot simultaneously diagonalize both matrices that correspond to the 
spatial differentiation along any two linearly independent directions. Therefore, to implement any 
characteristic-based treatment one has to select some specific direction for the spatial differentiation, 
which is equivalent to assuming that the external solution is composed of only some particular kind 
of waves. This immediately makes the characteristic-based treatment only approximate. 

The example of development and implementation of characteristic radiation boundary condi- 
tions for multidimensional problems is given in work by Thompson [63]. He studies quasi-one- 
dimensional formulation of the equations locally at each boundary point; the direction for the 
spatial differentiation is every time chosen to be normal to the boundary. This, in particular, 
implies that the exterior solution is assumed to consist only of the waves that propagate outward 
normally to the boundary. The ABC’s of [63] are constructed using the same type of character- 
istic analysis as described in [62] but applied here to the quasi-one-dimensional problem. In so 
doing, one cannot generally expect the boundary to be transparent for the outgoing waves with 
incident angles that differ from t /2. The approach of [63] was further developed and generalized 
by Vanajakshi, Thompson, and Black [64] and Thompson [65]. 

Generally, the approach based on local characteristic analysis is attractive for practical com- 
puting because of its algorithmic simplicity, low computational cost, and geometric universality. 
The last feature in multidimensional problems is achieved by introducing the foregoing essential 
simplification: the pre-selection of one specific direction for the spatial differentiation. This, at 
the same time, is a serious drawback of the characteristic-based approach, which eventually leads 
to relatively low final accuracy. In other words, the outgoing waves that propagate at an angle 
to the boundary and should therefore have been taken into account for specifying the incoming 
values at later times are partially reflected back to the domain, which effectively causes the wrong 
data coming from the boundary, and partially get through but never used then, which causes the 
irreversible loss of the corresponding information. Moreover, as mentioned by Kreiss an Gustafsson 
in [66] (see also [5]), the simplest radiation conditions do not apply to the case of nonzero forcing 
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terms/initial data outside the computational domain. Some other examples of flows to which radi- 
ation conditions do not apply because the incoming waves do exist in the solution can be found in 
the work by Thompson [63]. 

Finally, it is interesting to note that the radiation boundary conditions for the linear case can 
also be obtained using some low-order approximations to the symbols of the corresponding vpDO’s. 

Basically, it appears that the Euler radiation conditions based on (quasi- )one-dimensional char- 
acteristic analysis may provide sufficient accuracy only if the artificial boundary is located far 
enough from the source of perturbations (e.g., from the immersed body for external flow prob- 
lems). This circumstance, of course, may require an excessively large computational domain in 
comparison with that required by exact ABC’s. For each specific case, it is always a separate 
question whether or not the extra computational effort originating from this extended domain can 
be compensated for by the simplicity and low cost of the local ABC algorithm. 

One can also formulate some local boundary conditions for the Navier- Stokes equations, which 
form a system of the so-called incompletely parabolic type. A thorough analysis of well-posedness 
for some linear incompletely parabolic problems (in particular, for the linearized Navier-Stokes 
equations and for the linearized shallow water equations) is provided by Gustafsson and Sundstrom 
[67]. This analysis uses the apparatus of energy estimates and allows the authors to select appropri- 
ate (i.e., well-posed) boundary conditions in a certain initially prescribed class, namely, in the class 
of first-order differential relations. Later, Nordstrom [68] provided the same type of analysis for a 
wider class of boundary conditions — second-order differential relations. Local ABC’s derived in 
[67] and [68] can be successfully used for flow computations under the same restriction as applies to 
the characteristic Eider conditions: the computational domain must be large enough. We also note 
that the Euler radiation boundary conditions can be obtained from the constructions of [67, 68] as 
the Reynolds number vanishes. 

In the context of viscous flow computations we should also mention earlier work of Rudy and 
Strikwerda [69, 70], in which a special “hyperbolic-type” local NRBC was derived and numerically 
implemented for calculating the flow over a flat plate in the framework of the two-dimensional 
compressible Navier-Stokes equations. The study by Rudy and Strikwerda has probably for the 
first time revealed the superiority of NRBC’s over the simplest Dirichlet-type outflow boundary 
conditions for the Navier-Stokes computations. 

In later work [71], Abarbanel, Bayliss, and Lustman used linearization of the Navier-Stokes 
equations against the approximate wake-type downstream solution to construct local ABC’s for 
the viscous flow over a flat plate. The ABC’s [71] are based on selecting the long-wave modes from 
asymptotic expansion of the solution to the corresponding linearized system. The approach of [71] 
was then generalized by Danowitz [72] for a wider class of external viscous flows than only the flows 
over a flat plate. 

Another technique for constructing local ABC’s for viscous flow computations is based on 
extrapolation of flow variables from inside the computational domain to the artificial boundary; 
this technique can probably be referred to among the simplest ones. Generally, one cannot expect 
to obtain high accuracy from such a procedure. Moreover, for subsonic flows extrapolation may 
even lead to the ill-posed problems, as was (experimentally) demonstrated by Rudy and Strikwerda 
in [70]. However, for a particular class of problems, namely for the flows of the boundary layer type, 
when the solution contains strong transversal gradients, the extrapolation boundary conditions may 
be appropriate. Gustafsson and Nordstrom [73] provide some theoretical justification for using 
the outflow extrapolation boundary conditions to calculate steady flows with strong transversal 
gradients. Moreover, they corroborate numerically that if the entire outflow boundary is subsonic, 
then the extrapolation conditions are inapplicable; however, if the external flow is supersonic and 
the subsonic part of the flow is contained only in the boundary layer (which implies presence of the 
strong transversal gradients), then the extrapolation conditions do apply and produce good results. 
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The approach of [73] was further generalized by Nordstrom [74] for the case of unsteady flows and 
[75] for the finite-difference formulation of the problem. Similar issues have been then analyzed by 
Nordstrom in [76]. We should note that the case of intersection of the artificial boundary and the 
solid boundary (which means that the artificial boundary crosses the boundary layer) presents one 
of the most difficult situations from the standpoint of constructing the ABC’s. In the meantime, 
we do not know of any other robust technique for setting the ABC’s in this case except for the 
one based on extrapolation. We also note that the basic function of the extrapolation boundary 
conditions [73] is to suppress boundary layers that are generally relevant to the near-boundary 
behavior of the solutions to the Navier-Stokes equations. The idea of suppressing the non-physical 
boundary layers that may arise near the inflow and outflow artificial boundaries is a central one 
for the work of Johansson [77], in which the set of ABC’s that contain higher-order extrapolation 
conditions as a part was derived and numerically implemented for some incompressible Navier- 
Stokes computations. 

The basic conclusion that can be drawn from this introductory survey is that any algorithm for 
setting the ABC’s should, generally speaking, satisfy two groups of requirements, which to a certain 
extent may be contradictory. The first group is mainly associated with accuracy requirements; as 
mentioned above, meeting these requirements virtually always leads to the nonlocal ABC’s. On 
the other hand, the exact nonlocal ABC’s have their essential limitations. Namely, such boundary 
conditions can be accurately derived only for the linear governing equations (in most cases) and 
for some particular geometries; moreover, the discrete implementation of nonlocal ABC’s is not 
always easy, even when their continuous formulation is available. The second group combines 
the requirements of algorithmic simplicity, geometric universality, and low computational cost; 
the approximate approaches that result in local ABC’s are usually much better than the exact 
ones from this standpoint. On the other hand, the accuracy provided by local ABC’s is often 
insufficient. Therefore, in constructing a specific numerical algorithm, one always has to select an 
optimal computational strategy compromising between the two groups of requirements mentioned 
above. 

However, modern trends in the development of numerical methods will, in our opinion, make 
higher and higher requirements for the accuracy of computations. Consequently, we expect the 
practical demands for highly accurate ABC’s to grow in the future. This will obviously necessitate 
paying more attention to constructing such highly accurate (nonlocal) ABC’s that will at the same 
time be computationally effective (a few particular examples of this type are described below). 
Additionally, one should always aim to obtaining a robust numerical procedure, therefore, the 
robustness becomes one of the most important experimental criteria for assessing and comparing 
different ABC’s. 

Obviously, the review provided above is not complete. Referring the reader to two comprehen- 
sive works by Givoli [78, 79] for more survey information, we now proceed to describing the DPM 
[29, 30] and the DPM-based ABC’s, which are the focus of this presentation. 

2. The Generalized Potentials and the Difference Potentials Method. Our ultimate 
goal is to construct ABC’s that, to a certain extent, combine the advantages of local and nonlocal 
approaches. Namely, we aim to achieve high accuracy relevant to the nonlocal techniques and, at 
the same time, geometric universality and algorithmic simplicity relevant to many local methods. 
Our main tool in achieving this goal is the apparatus of the generalized potentials and the DPM 
[29, 30]. 

In this section, we use a model example for the Poisson equation to show the principle elements 
of the construction of generalized potentials and their implementation for setting the ABC’s. Of 
course, the presentation below is far from being thorough. A delineation of the ideas associated with 
generalized potentials, boundary equations with projections, and their numerical implementation 
(DPM), as well as many examples, can be found in the original work by Ryaben’kii [29, 30]. 
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The next section of this paper will be devoted to the specific constructions of ABC’s based on 
the DPM. 

Let us first return for a moment to the second example in Section L There, we have truncated 
the original infinite problem (1.3) and have obtained the finite problem (1.3a)-(1.5) on the disk 
{r < i?o} instead. The artificial boundary {r = i2 0 } has been chosen to be circular, which is 
essential for constructing ABC’s (1.5) (indeed, the separation of variables could not be implemented 
otherwise). 

Here, we are going to consider an analogous example but for three space dimensions. Switching 
from two to three dimensions preserves all essential elements of our construction and, at the same 
time, allows us to avoid the unnecessary complications associated with the logarithmic asymptotics 
of some classical potentials in two dimensions. Thus, let us consider on R 3 the Poisson equation 
with the compactly supported RHS: 


(2.1a) An = f(x) y x € R 3 , supp/(x) C B C {\x\ < R 0 }. 

We impose the condition of vanishing of the solution to (2.1a) at infinity, 


(2.1b) u(x) — * 0, as |x| — * oo, 

and note that unlike the 2D case, we do not need to demand any special properties of f(x) (compare 
with (1.3a)) to ensure the solvability of (2.1). Expanding the solution and the RHS in terms of 
spherical functions, we obtain for |x| > Ro the following family of ordinary differential equations 
describing the radial modes (r = |x|): 


(2.2) i.( r 2^l) _/(/ + 1)^ = 0 , / = 0,1,2,-... 

As in the 2D case, equation (2.2) has two eigensolutions for each /, l = 0,1,2,...: «^(r) = r l 
(increasing) and t*| 2 ^(r) = (decreasing). To satisfy (2.1b), we need to select only the 

vanishing mode from the above two. Therefore, the following set of relations (compare with (1.5)) 


(2.3) 


diii 

dr 


i ^ + * 

H u\ 

t—Rq t 


r=Ro 


— 0 , / = 0 , 1 , 2 , . 


provides here the exact ABC’s at the spherical artificial boundary {r = Ro}- (Note, 21 + 1 linearly 
independent spherical functions correspond to each value of /; this circumstance, in a sense, does 
not influence the construction of ABC’s (2.3) since the boundary conditions are simply the same 
for all 21 + 1 components.) 

The problem obviously becomes much more difficult if, for some reason, we need to consider an 
artificial boundary with a more complicated shape. In practice, the shape of an artificial boundary 
is often prescribed by the algorithm (more precisely, by the grid) used inside the computational 
domain; some real examples will be provided in Section 3. As concerns the model case studied 
here, we simply assume that there is an irregular artificial boundary T that separates the finite 
computational domain D{ n from its infinite exterior D ex ; we also assume that Di n entirely contains 
the support B of the RHS f(x) (see (2.1a)). The geometric setup (projection onto the plane) 
is schematically shown in Figure 2.1, in which the closed dashed line represents T. The sphere 
{ r = R 0 } is needed as an artificial element in further consideration; without loss of generality, we 
may always assume that D{ n C {r < i?o}- 
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FlG. 2.1. Example for Poisson equation in three dimensions: schematic geometric setup (projection onto the plane). 

We now recollect the classical Green formula for harmonic functions. Namely, let u ex (x ), 
x G R 3 y be harmonic in D ex and have a zero limit at infinity, u ex (x) — * 0 as |x| — » oo. Then, 


(2.4) « ex (x) = -©r(x)- 


c(*) 


+ 


x( 


dG(x,y) 


Uex(y) “ G{x,y) 


du ex (y) 


ds 




x e D e 


where ©r(x) is the characteristic function of set T, ©r(^) is equal to one on T and zero on R 3 \T; 
G(x,y) — — (dx) -1 ^ - y I” 1 ; n is the outward normal to T; s is the surface area; and the subscript 
y denotes the variable of differentiation or integration. Note, the first term on the right-hand side 
of (2.4) takes into account the jump of the double-layer potential, the addition of this term makes 
(2.4) valid on the closed domain D ex . 

We emphasize that representation of u tx (x) as a sum of the two potentials (double-layer and 
single-layer, see (2.4)) holds only for harmonic functions on the domain. If, however, we specify 
two arbitrary functions on T and substitute them into (2.4) as densities of the potentials, then the 
resulting function will obviously be harmonic on jD ex , but its boundary values, as well as boundary 
values of its normal derivative, will not, generally speaking, coincide with the original densities of 
the double-layer and the single-layer potentials, respectively. 

Let us now define the generalized potential with vector density £r = (fo?fi) specified on T. 
We will use the formula analogous to (2.4) but will not require in advance that £ 0 and be the 
boundary values of some harmonic function and its normal derivative, respectively. Specifically, 
the generalized potential is given by the following expression: 


(2.5) u ex (x) = P ei£r = -Qr(j) ^°i - + 


/( 


dGj^y) 

dn v 


£o (y) - G(x,y)^(y) ds y , x 6 D e 


We also define the operation of taking the (vector) boundary trace of the function defined on D ex : 
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(2.6) fr = Tr u ex (x) = (u ex , ^^) . 

Finally, we define the boundary operator Pp that maps the space of traces £ r onto itself; this 
operator is a composition of the generalized potential P ex and the trace Tr, 


(2.7) P r fr = TVP c *fr. 

Referring the reader to [29, 30] for further details, we only note here that the structure of the 
operator Pp from (2.7) is nontrivial since, in particular, it contains the normal derivative of the 
double-layer potential. The corresponding singularity may noticeably complicate the direct calcu- 
lation of Pp£p; therefore, the actual calculation requires the development of a special alternative 
approach (see below). However, the operator Pp itself appears to be extremely useful in the analysis 
of boundary- value problems, it plays a fundamental role in our further consideration. 

We note that the operators analogous to Pp for the first time appear in the work of Calderon 
[80] and then Seeley [81]. It turns out that Pp is a projection, Pr = P^; sometimes it is called the 
Calderon boundary projection. Ryaben’kii (see [29, 30] and the corresponding bibliographies) had 
modified and generalized the original construction and developed an effective way for numerical 
treatment of the boundary projection operators. 

It is easy to show that those and only those vector functions £p that satisfy the following 
boundary equation with projection (BEP) 


( 2 . 8 ) 


Pr£r = fr 


admit a harmonic complement u ex on D ex that vanishes at infinity and has the trace fp on T, 
Tr u tx = fp. Indeed, if u ex (x) is harmonic on D ex and vanishes at infinity, then we rewrite Green’s 

OV/QX 


formula (2.4) as u ex (x) = P e 


dn 


and, applying Tr from (2.6), we obtain (2.8). Con- 


versely, if (2.8) holds for some fp, then the harmonic function P ex fp vanishes at infinity (as a sum 
of two classical potentials) and has the trace £p. 

Thus, equation (2.8) is paramount since it provides an exhaustive classification of those and only 
those vector densities £p that have a harmonic continuation on D ex . In other words, equation (2.8) 
is equivalent to the Laplace equation on D ex along with the condition of vanishing of its solution at 
infinity. Therefore, equation (2.8) provides us with a desirable exact ABC at the irregular boundary 
T. Rather than solving (2.1) or (2.1a)-(2.3), we may now solve (2.1a)-(2.8) and obtain exactly the 
same solution on D { n . 

The operator Pp is obviously nonlocal. Moreover, some computationally difficult elements still 
remain in its structure. Therefore, we will now redefine Pp in a more practical way. First, let us 
recall the classical Green identity (a is the volume): 


(2.9) 


-0r(z) — - — + 


dG(x,y ) , ,du ex (y) 

SZ u cx(y) ~ G(x, y)— — 


dSy — 


= u ex (x ) - J J D G{x,y)AyU ex (y)da y , x € D ex , 


which holds for any function u ex (x) such that the convolution in the second term on the right-hand 
side of (2.9) exists. Since (2.9) is an identity, and its left-hand side coincides with the Green formula 
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for harmonic functions (see (2.4)), we can use (2.9) to construct a new definition of generalized 
potential. Indeed, for any prescribed vector density £r> we simply take some u(y) defined on R 3 
and such that Tru = Then, we consider the following function: 


( 2 . 10 ) 


9ex(y) — 


Ayu(y), y e D ex 
o, ye Din 


and easily derive 


(2.11) P e *£r = u( x ) - J J R3 G{x, y)g ex (y)da y , x € D ex . 

Indeed, formula (2.11) immediately follows from definition (2.5) and identity (2.9) since, according 
to (2.10), the right-hand sides in (2.9) and (2.11) are the same for x € D ex . As concerns the choice 
of u(y ) from (2.10), the only essential requirement is that Tru = £p; this requirement is always 
easy to meet. At the same time, we can easily ensure the existence of convolution from (2.11), since 
the function u(y) (and, consequently, g e x(y)) can, for example, always be chosen to be compactly 
supported without loss of generality. 

Using this new definition of generalized potential (see (2.11)), we construct the boundary 
projection Pr in the same way as before, i.e., in accordance with formula (2.7). We emphasize that 
the new definitions of P ex and Pp no longer require the calculation of surface integrals; instead, 
one needs to calculate a volume Newton potential (see (2.11)). 

The last step is to understand that, generally speaking, we do not need to know the generalized 
potential everywhere on the infinite domain D ex . In particular, to construct exact ABC’s in the 
form (2.8) we need to know this potential only in the vicinity of T to calculate the operation Pp. 
Therefore, we can, for example, introduce the spherical artificial boundary {r = Rq} and rewrite 

(2.11) as 


(2.12) Pexfr = M>) - / / G(x, y)g ex (y)d<T y , x € D ex n {r < R 0 }, 

J J{t<Rq} 

assuming (without loss of generality) that supp g ex (y) C {r < iio}* (Recall, u{y) from (2.10) can 
always be chosen to be compactly supported.) Clearly, boundary projections (operators Pp) that 
can be constructed on the basis of (2.11) and (2.12) are exactly the same. Now, we simply notice 
that the second term on the right-hand side of (2.12) is the solution of a certain auxiliary problem 
(AP) formulated on the ball {r < i?o}* Indeed, this term obviously solves on {r < R 0 } the Poisson 
equation 


(2.13) jAu — ^ex(*^)> ^ ? 

supplemented by boundary conditions (2.3). Therefore, AP (2.13)-(2.3) can be used instead of 
(2.12) for calculating the generalized potential on {r < # 0 } for this specific example. Usage of 
the AP (2.13)-(2.3) instead of convolution (2.12) is a major simplification from an algorithmic 
standpoint since the AP is formulated on a simple domain and admits an effective solution by the 
separation of variables. 

Generally, the concept of the AP is the second most important element in our consideration. 
Indeed, we see that exact ABC’s at the irregular artificial boundary T are obtained in the form of 
BEP (2.8). To calculate the projection Pp, we originally needed surface integrals (2.5), then we 
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replaced the surface integrals with a volume potential (2.11), and finally, we have found the way to 
calculate the generalized potential and the boundary projection using AP (2.13)-(2.3), which can 
be easily solved by separation of variables. 

Summarizing, we see that to set the exact ABC’s at the irregular artificial boundary one needs 
first to formulate an appropriate AP and then to construct the BEP of type (2.8). We do not further 
delineate this statement here; instead, we provide some practical examples in the next section. 

To complete this section, we first note that metric properties of potentials and projections 
(well-posedness of BEP’s) for different problems are thoroughly studied by Ryaben’kii in [29, 30]. 
Second, we briefly address the numerical aspects of calculating generalized potentials and projec- 
tions, i.e., the DPM in a narrow sense. The essence of the DPM is an analogue of the Calderon 
boundary projection constructed for the finite- difference formulation. We begin with specifying 
some difference AP; the AP should be uniquely solvable and well-posed and should also admit an 
easy numerical treatment. For example, AP’s that can be solved by separation of variables on the 
basis of, e.g., the discrete Fourier transform or the difference spherical functions (see [7]), have 
found extensive practical applications (see [29] and the examples below). Of course, the difference 
AP used for constructing the ABC’s should contain some element analogous to (1.5) or (2.3) for the 
exact transfer of boundary conditions from infinity, examples of such AP’s are given in Section 3. 
The grid for the AP is usually regular (e.g., Cartesian, polar, spherical, etc.), and no grid adaptation 
to the shape of T is required. Instead, we consider a special subset of nodes of this grid called the 
grid boundary. The grid boundary 7 is located, in a certain sense, close to the continuous boundary 
T; the structure of 7 depends on the stencil used for the difference approximation of the original 
differential operator. Specifically, the entire regular grid is separated by T into two subsets of nodes; 
the first subset belongs to D tn , and the second one belongs to D ex . These subsets have no common 
points. We now apply the stencil to every node of each of the two above-mentioned subsets; clearly, 
the stencil always sweeps a wider grid area than the original subset. The two grid areas swept by the 
stencil obviously have now a non-empty intersection; the latter is called grid boundary 7. Basically, 
it appears that 7 is a multi-layered set of nodes of the regular grid concentrated near T. Instead 
of vector densities £r? we consider scalar grid functions on 7; since 7 is multi-layered, the former 
can in a certain sense be modeled by the latter. The formal construction of generalized difference 
potentials and difference projections consists of the same elements as above. We simply plug the 
corresponding difference objects instead of the continuous ones into the aforementioned procedure. 
For example, the corresponding subsets of the regular grid are substituted for D in and D tx \ the grid 
boundary 7, for T; the difference operators, for the differential ones; solution to the difference AP, 
for the solution to the continuous AP. In so doing, we obtain difference analogues to the potentials 
and projections; the issues of consistency and convergence are delineated by Ryaben’kii in [29]. 
The full scheme of the DPM (and the construction of the DPM-based ABC’s) may also require the 
operation of continuing the boundary data from T to 7, as well as some interpolation operations. 
Continuation of data from T to 7 is usually done by means of the Taylor formula, the details can 
be found in [29, 30]. Difference potentials and projections, which can be effectively computed in 
practice, are used for solving boundary- value problems by means of the DPM, as well as for setting 
the ABC’s. 

Finally, we note that implementation of the DPM for setting the ABC’s should not be assessed 
only from the standpoint of geometric universality. It comes out that only some particular classes of 
PDE’s admit a true exact transfer of boundary conditions from infinity to a closed finite boundary 
(e.g., circle or sphere). Indeed, such a transfer generally requires some sort of symmetry (e.g., 
cylindrical or spherical), which holds, for example, for the wave, Laplace, or Helmholtz equation, 
but cannot be referred to as a general case. For the general case, even the construction of a 
closed finite boundary T, which would enable us to exactly transfer the boundary conditions from 
infinity (or to do that as close to exact as desired), may require some sophisticated elements (see 
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below). In our opinion, such constructions, which are naturally incorporated in formulations of the 
corresponding AP’s, can be most effectively implemented in just the framework of the DPM. 

3. The DPM-Based Nonlocal ABC’s. General theoretical foundations for constructing 
the exact DPM-based ABC’s are provided by Ryaben’kii [17], see also brief note [82], This work 
presents an important contribution to the literature on numerical methods for problems on un- 
bounded domains since the BEP’s of type (2.8) are for the first time used in [17, 82] as the ABC’s 
to be set when truncating an external problem for the purpose of its numerical solution. Formerly, 
the Calderon boundary projections were mostly applied to the qualitative analysis of integral equa- 
tions and PDE’s, and the DPM was mostly implemented only for solving the boundary-value 
problems. 

The author of [17, 82] studies a general unsteady problem formulated on some infinite domain. 
The problem is assumed to be already discretized on some grid, so that only the finite-difference 
formulation is considered. This allows one to obtain the ABC’s directly for the specific numerical 
algorithm rather than for the original continuous formulation, which is convenient in practice. 
Outside some finite grid domain (i.e. , the computational domain), the problem is assumed to 
be linear and homogeneous, and the original boundary conditions at infinity are also assumed 
to be homogeneous. No restrictive assumptions are made in regard to the problem inside the 
computational domain; it can even be nonlinear. The only important requirement of the overall 
formulation is that the problem be uniquely solvable and well-posed. Such a consideration is relevant 
to many applications, for example, solid mechanics, where one often has some bounded region with 
the strong deformations (like plasticity or destruction) surrounded by an extended medium where 
the deformations are small and are, therefore, governed by the linear elasticity equations. In regard 
to the shape of the computational domain in [17], generally speaking, no restrictive requirements 
re imposed. 

The ABC’s [17] are obtained in the form of a difference BEP analogous to (2.8). As mentioned 
above, the difference BEP’s are written with respect to a grid potential density defined on the grid 
boundary 7 . In the case of unsteady problems, 7 is generally a multi-layered cylindrical surface (with 
an element aligned to the time axis) consisting of grid nodes; it is an analogue to the continuous 
cylindrical boundary of the P X [0, T] type. The AP used by Ryaben’kii in [17] is also unsteady; 
it is solved by some time-evolution technique. Since the Green operator of AP is incorporated 
in the structure of the boundary projection, the resulting ABC’s in [17] appear to be nonlocal in 
both space and time. We emphasize that these ABC’s are exact in the sense mentioned above: the 
solution found inside the computational domain with the help of these boundary conditions is the 
same as it would be if the original (infinite-domain) problem were solved. 

To reduce the computational cost associated with nonlocality, which can be rather high for 
the general formulation, Ryaben’kii [17] proposes several different approaches that are relevant to 
some particular classes of problems. For example, if the coefficients of the linear system, as well 
as the discretization parameters in space, do not depend on time, then the corresponding Green 
operator becomes invariant with respect to a shift along the time axis. This provides a noticeable 
economy when calculating and storing the coefficients of the Green operator. If the problem under 
study is parabolic, then the coefficients of its Green operator become small as the corresponding 
time interval becomes large. Therefore, only some fixed number of these coefficients is effectively 
needed to be taken into account for each specific moment of time, whereas the remainder of the 
coefficients can be neglected. This obviously leaves the ABC’s nonlocal in time but only for some 
finite interval in the past. For some hyperbolic problems, it is also possible to effectively “cut off 
the tail” in time using the existence of lacunas. In particular, if the number of space dimensions 
is odd, the coefficients of a linear hyperbolic system are constant, and the lower order terms in 
the equations are absent, then the solution to the corresponding Cauchy problem has lacunas. In 
particular, it means that if the initial data are compactly supported, then the solution becomes zero 
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in any fixed (spatial) domain for sufficiently large time. The same property should asymptotically 
hold for the difference Cauchy problem as well; the latter can be used as an AP. In so doing, we 
again can neglect those coefficients of the Green operator that correspond to the large time intervals 
and, therefore, restrict the nonlocality of the ABC’s in time. These techniques, as well as some 
others presented in [17], can substantially reduce the required computer effort when calculating 
the nonlocal exact ABC’s. Finally, it is worth mentioning that the DPM-based ABC’s [17] can be 
considered as a certain type of “standard software”. Indeed, the parameters of the process inside 
the computational domain can be changed, but as long as the exterior linear problem and the shape 
of the artificial boundary remain intact, the ABC’s also remain exactly the same. Therefore, these 
boundary conditions can be effectively used for solving families of similar problems. 

It is important to note that from an algorithmic standpoint, the role of the ABC’s (in particular, 
of type [17]) is simply to close the finite- difference system solved inside the computational domain. 
Indeed, the stencil of the scheme used inside D{ n cannot be applied to the nodes located at the 
external boundary or near it since the part of the stencil may simply “fall out” of the domain. 
Therefore, the finite- difference system appears to be undetermined and requires that the missing 
relations between the values of the solution near the external boundary be provided by the ABC’s. 
When the ABC’s are constructed as a difference analogue to (2.8), one can always think that the 
data on the internal layer(s) of the grid boundary 7 are known and the data on the external layer(s) 
are to be determined using the BEP. Therefore, the question of actually solving the BEP (for a 
specific finite-difference formulation) and expressing the boundary values in terms of data provided 
from inside D xn becomes an important issue. This question is closely related to another one, namely, 
implementation of the DPM-based ABC’s with different types of time-evolution procedures. Some 
examples of such an implementation are also given in [17], 

Ryaben’kii [17] provides a general framework of the algorithm for setting the ABC’s based 
on the DPM. Any specific problem, of course, requires special treatment. Zueva, Mikhailova, and 
Ryaben’kii [83] and Brushlinskii, Ryaben’kii, and Tuzova [84] numerically investigate diffusion of 
the magnetic field in a finite conducting object immersed in a vacuum. Both papers deal with the 
2D formulation. In [83], the Cartesian coordinates are used to study the magnetic field on the 
plane normal to a finite-section conducting rod; in [84], the conducting structure is a finite- section 
axisymmetric annulus, and the magnetic field is studied in the meridian plane (r, z). The diffusion 
process inside the conducting structure may be unsteady; however, the field in the surrounding 
vacuum is described by the Laplace equation. Therefore, the AP in [83] and [84] is formulated 
and solved for the Laplace equation. Physical formulation of the problem generally assumes that 
the surrounding vacuum is infinite; in practice, however, the infinity is modeled by some remote 
boundary at which the magnetic field turns to zero. Clearly, the size of the entire domain (including 
the vacuum area) must be chosen to be much larger than the cross section of the aforementioned 
rod or annulus; consequently, the ABC’s that can exactly transfer the remote boundary condition 
to the surface of the conducting structure are as important in this model formulation as in the truly 
infinite one. In [83, 84], the remote boundary condition is simply incorporated into the formulation 
of the AP, Then, the exact ABC’s are obtained in the form of a difference BEP; the BEP is further 
resolved for numerical convenience. This process results in a matrix relation that connects the 
vectors of unknowns on internal and external layers of the two-layered grid boundary 7 . Therefore, 
the implementation of the ABC’s [83, 84] becomes simple from an algorithmic standpoint: the 
matrix is calculated only once, and then it multiplies the corresponding vector at each time step. 

The DPM-based ABC’s [83, 84] can be obtained in exactly the same manner for any shape 
of the cross section of the conducting object; the actual numerical experiments in [83, 84] are 
carried out for the rectangular cross section. These experiments show that the computational cost 
of the algorithm that includes the DPM-based ABC’s is much lower than the cost of the standard 
procedure that is based on actually solving the Laplace equation in a vacuum area at each time 
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step. At the same time, the accuracy of the solution calculated inside the rod (annulus) with the 
ABC’s [83, 84] is the same as the accuracy obtained when the entire original problem is solved on 
an extended domain. Therefore, the DPM-based algorithm turns out to be an attractive technique 
for solving problems of this type. 

In [85], Mishkov and Ryaben’kii solve the nonhomogeneous Helmholtz equation on a semi- 
plane. The RHS of the equation is compactly supported, it models a sound source (e.g., a scatterer 
or an emitter). The surrounding semi-infinite sound-conducting medium is stratified; it consists 
of two layers with different speeds of sound, which makes the coefficient at the zeroth-order term 
piecewise constant. Finally, radiation boundary conditions at infinity are formulated using the 
limitary absorption principle. The solution to this problem should be numerically found on some 
finite computational domain that entirely contains the aforementioned sound source. 

The non-constancy of coefficients makes this problem difficult to handle. In particular, there is 
most likely no means to exactly (analytically) transfer the boundary conditions from infinity even to 
some regular finite boundary, as in the case of the Laplace equation that was done for a circle and for 
a sphere (see (1.5) and (2.3), respectively). Therefore, the AP in [85] is formulated for the perturbed 
Helmholtz equation (finite absorption is added) on a sufficiently large rectangle, and zero Dirichlet 
boundary conditions are set at its sides. It is shown that as the size of this rectangle enlarges, 
the solution to the AP approaches the solution of the corresponding infinite-domain problem on 
any fixed neighborhood of the computational domain. It is also shown that as the absorption 
coefficient vanishes, the solution to the difference AP does approach the solution that corresponds 
to the true outgoing waves without an absorption. Moreover, a special method is proposed in [85] 
that enables one to formally let the absorption coefficient be zero in the computations and to still 
obtain the solution composed of only outgoing waves. Consequently, ABC’s [85] that are obtained 
in the form of a difference BEP based on the solution to the difference AP described above, can 
generally be constructed as close to the exact ABC’s as desired. In other words, with these ABC’s 
the truncated problem can be solved so that its solution differs from the corresponding fragment of 
the original solution within the initially prescribed accuracy. Of course, achieving the high accuracy 
may require a large domain for the AP. However, ABC’s [85] become most effective when they are 
used to solve a series of similar problems, e.g., to calculate sound fields from different sources in 
the same medium, which is frequently an important practical case. 

Tsynkov [86, 87] and Sofronov and Tsynkov [88] consider a finite body (airfoil) immersed 
in an infinite flow of inviscid compressible fluid. The flow is governed by the Euler equations 
and is assumed to be subsonic at infinity; for the purpose of numerical solution, the equations are 
discretized on a finite-difference O-type grid that is generated around the body. The computational 
domain in [86, 87, 88] is formed by the grid; therefore, the shape of its external boundary is 
completely determined by the grid as well, and no special assumptions in regard to this shape 
are made. Outside the computational domain, the Euler equations are linearized against the free- 
stream background. Moreover, we assume the existence of the velocity potential in the far field 
and split the linearized system into elliptic (velocity) and advection (entropy) parts. After the 
term associated with the circulation of flow around the airfoil is subtracted, the regular part of the 
potential of velocity perturbations is described by the Prandtl-Glauert equation 

( 31a ) ( l ~ M ° 2 ) + = 0 
with a zero boundary condition at infinity: 

(3.1b) <p{x,y) — ► o, as x 2 + y 2 — ► oo. 

Equation (3.1a) can be easily reduced to the Laplace equation by means of an affine coordinate 
transformation; in so doing, boundary condition (3.1b) obviously remains intact. To obtain ABC’s 
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for velocity components, we represent the solution to (3.1) in the form of a generalized potential 
and then construct the corresponding BEP. We use the AP formulated on a ring-shaped domain 
{Ri < r < i? 0 }; the external circle {r = R 0 } encompasses the artificial boundary and the internal 
circle {r = R x } lies inside the computational domain. At r = R u we specify homogeneous Dirichlet 
boundary conditions; at r = R 0 , we specify boundary conditions (1.5a). It is possible to make 
sure that such an AP is uniquely solvable and well-posed for any RHS concentrated inside the 
ring. Numerically, the AP is easy to solve by means of the discrete Fourier transform in polar 
coordinates. 

As mentioned above, the ABC’s in the discrete formulation should simply close the system that 
is solved inside the computational domain. For the case of a second-order scheme employed inside 
the computational domain, we can always assume that velocity components on the penultimate 
coordinate row T of the O-type grid are known, whereas the corresponding values on the outermost 
row r i should be determined using the ABC’s. Therefore, we consider the penultimate coordinate 
line T as an actual artificial boundary. Grid boundary 7 consists of those nodes of the polar AP 
grid that are located near T (see Section 2). The finite-difference BEP with respect to a potential 
density defined on 7 is constructed in a standard way (see Section 2) on the basis of the AP 
described above. Then, this BEP is solved so that the velocity on T provided from inside the 
computational domain coincides with the gradient of potential (3.1) in a certain generalized sense. 
After the BEP is solved, we can find the grid potential density for any data (velocity components) 
specified on T. Once this grid density is known, we calculate the generalized potential and find the 
trace of its gradient on the outermost coordinate line Ti by means of interpolation; this procedure 
yields the ABC’s for velocities. Finally, the ABC’s for thermodynamic parameters are obtained 
using local relations, specifically, the Bernoulli equation and the entropy advection equation. 

The technique of [86, 87, 88] for constructing the ABC’s was combined with an iterative method 
[89] by Sofronov for calculating steady solutions to the Euler equations. A few subsonic and 
transonic airfoil flows have been numerically studied; in the transonic case, local supercritical 
regions should always be located inside the computational domain so that the exterior linearized flow 
remains purely subsonic. Results of the numerical experiments presented in [88] demonstrate clear 
superiority of the nonlocal DPM-based ABC’s over the standard local techniques based on quasi- 
one- dimensional characteristic analysis (see above). For a fixed computational domain, nonlocal 
ABC’s [86, 87, 88] provide better accuracy and a faster convergence rate than local techniques; the 
entire numerical algorithm also appears to be more robust. Indeed, for some variants the iteration 
procedure supplemented by local boundary conditions may simply fail to converge, which never 
occurs for ABC’s [86, 87, 88]. Additionally, when the artificial boundary approaches the airfoil, the 
solution obtained with the technique of [86, 87, 88] appears to be essentially less influenced by the 
decrease in the size of the computational domain than the solution obtained on the basis of the 
local boundary conditions. In other words, ABC’s [86, 87, 88] allow one to maintain good accuracy 
for much smaller computational domains than the standard boundary conditions do. These results, 
along with the geometric universality of ABC’s [86, 87, 88], make this approach most useful for 
calculating external Euler flows. 

Let us now describe another algorithm for setting the nonlocal DPM-based ABC’s in external 
flow computations. The foundations of this algorithm are proposed by Ryaben’kii and Tsynkov 
[25], further generalizations, as well as results of the numerical experiments, can be found in the 
work of Tsynkov [1] and Tsynkov, Turkel, and Abarbanel [90]. 

We consider a finite body immersed in an infinite flow of viscous compressible fluid; at infinity, 
the flow is assumed to be subsonic. The general framework for constructing the ABC’s remains 
the same as above. We first linearize the governing equations in the far field, then we formulate 
an appropriate AP, and finally we obtain the corresponding BEP. However, the treatment of the 
flow in a viscous formulation requires the development of special approaches for all stages of the 
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aforementioned procedure. Below, we concentrate primarily on the construction of the AP since 
the AP is responsible for the “right” far-field behavior of the solution; therefore, the AP is the main 
element of the entire ABC algorithm that reflects the essential features of the original formulation. 



Fig. 3.1. Configuration of domains. 


The general geometric setup for the problem is shown in Figure 3.1. We integrate the Navier- 
Stokes equations on a finite- difference C-type grid generated around the body (airfoil). This grid 
covers the finite computational domain D in , the infinite exterior to Z?, n is called D ex . The domains 
Di n and D ex are separated by the artificial boundary P, which is the penultimate coordinate row 
of the grid. The outermost coordinate row of the grid is called Tx. As mentioned above, the ABC’s 
should close the finite-difference system solved inside D tn . In other words, for a second-order scheme 
employed inside the computational domain (3x3 stencil) the ABC’s should provide the missing 
relations between the values of the solution at P and at Ti. We will further discuss this case only; 
if, however, the stencil used inside D{ n is larger than 3x3, then the entire construction requires 
only minor changes. Namely, instead of Pi one will have to consider more coordinate lines exterior 
to r. 

Assuming that the flow perturbations are small in the far field, we linearize the governing equa- 
tions in D ex around the constant free-stream background. The linearized stationary compressible 
Navier-Stokes equations (2D plane dimensionless formulation) take the form 

dp du dv _ 
dx dx + dy 

du dp 1 / 4 d 2 u 1 d 2 v # 2 tt\ _ 

dx dx Re y3d:r 2 + Zdxdy ^ dy 2 ) ’ 

dv dp 1 (a d 2 v 1 d 2 u d 2 v\ _ 

dx dy Re y3dt/ 2 + 3 dxdy dx 2 J 


(3.2a) 
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dp 1 dp k ( 1 \ 

dx ~ M 0 2 dy ~ OFF \ Ap ~ kMq 2 ^) ~ °’ 

where u, v, p, and p are the perturbations of velocities, pressure, and density with respect to the 
corresponding free-stream parameters; Mo is the Mach number at infinity; Re is the Reynolds 
number at infinity; Pr is the Prandtl number; and k is the ratio of specific heats. We use the 
equation of state (k - l)e = p/p to eliminate the internal energy e from the fourth equation in 
(3.2a); A in this equation is the Laplace operator. 

System (3.2a) is supplemented by the condition of vanishing of all perturbations at infinity: 


(3.2b) (u, v, p, p ) — » (0, 0, 0, 0), as x 2 + y 2 — ► oc. 

Further, instead of the original formulation we consider a new coupled problem; this problem is 
nonlinear in D tn (original Navier-Stokes equations with no-slip conditions at the body) and linear 
in D ex (system (3.2a) with the free-stream limit of the solution (3.2b) at infinity). We cannot 
solve this coupled problem directly because D ex is infinite, so we equivalently replace the entire 
linear part of the problem by the BEP formulated on T (more precisely, by the difference BEP). 
To accomplish this step, we need to formulate an AP, which takes into account the structure of the 
solution from outside the computational domain. Recall, in Section 2 such an AP was formulated 
on a ball and the boundary conditions were transferred from infinity to the surface of this ball using 
an analytical approach. For the example in Section 2 this approach w as feasible since the Laplace 
equation admitted the separation of variables in spherical coordinates. 

Here, we would like to develop an analogous construction for (3.2). However, most likely no 
analytical approach can be implemented directly to exactly transfer boundary condition (3.2b) 
for system (3.2a) from infinity to some closed finite boundary. Apparently, system (3.2a) admits 
separation of variables in the Cartesian coordinates only, which is a major difference in comparison 
with the cases investigated previously. Obviously, this circumstance will strongly influence the 
formulation of the AP. 

First, we introduce a rectangular domain Dy = (0, X) x (— F/2, Y/2) D (see Figure 3.1) 
and construct a central- difference second-order approximation to (3.2a) in Dy on a Cartesian grid 
with the sizes h x and h y . In accordance with the general scheme of the DPM, AP is formulated 
for the inhomogeneous difference counterpart to (3.2a). A compactly supported RHS for the AP is 
specified in the same way as described Section 2 (see [25] and [1] for more details). 

Further, we assume periodicity with the period Y along the y direction and take the discrete 
Fourier transform (with respect to y) of both parts of the inhomogeneous difference analogue to 
(3.2a). This process results in a second-order system of ordinary difference equations for each 
wavenumber k. Then, we reduce this system to a first-order one by introducing the new variables 
[25] and obtain 

(3-3) A*v° + B*v° _ u = g° 

m — 1, . . . , M, k — — J, 

where ~ denotes the Fourier transform and M + 1 and 2 J + 1 are the number of nodes of the 
Cartesian grid along x and y, respectively. Explicit expressions for the entries of the 8x8 matrices 
A* and B^, as well as the definitions of the eight-component vectors of unknowns k and RHS’s 
Sm k are gi yen m [25]. Boundary conditions for the AP at x = 0 and x = X are imposed on Fourier 
components of the solution separately for each wavenumber k as 
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(3.4a) II (Q*-M*(*)I)v§ ifc = 0, k = 0,±l,±2,...±J, 

(3-4b) El (Qfc-^(fc)I)vS, ifc = 0, * = 0,±1,±2,...±J, 

where Q* = AJ x Bfc in (3.4), I is the identity matrix and ^ s (fc), s = 1,. . .,8, are the eigenvalues 
of Qjt (to be found numerically in practice). 

Boundary conditions (3.4) are the principal part of the formulation of AP. Since the RHS is 
compactly supported, we may formally consider (3.3) for m < 0 and for m > M as a homogeneous 
system of ordinary difference equations, which has (exponentially) increasing, as well as decreasing, 
eigensolutions. Then, condition (3.4a) prohibits at x = 0 (m = 0) all modes that do not decrease 
to the left (inflow direction), and condition (3.4b) prohibits at x = X (m = M ) all modes that 
increase to the right (outflow direction), see [25]. 

We now introduce the following concept of convergence, (i) The continuous AP can be for- 
mulated in a natural manner analogous to the difference AP (see [25]). (ii) The convergence of 
a difference solution to a continuous one is considered only on any finite subset (rectangle) of the 
type (0, X ) X (— y, y) D D in , y < Yj 2, rather than on the whole Dy. (iii) We consider con- 
vergence not only while the grid size vanishes but also while the period Y synchronously grows 
(see [25]), namely, (h r , h y , Y) — > (0, 0, -foe). Note that the growth of the period Y implies 
that we proceed from periodic functions to nonperiodic ones; the latter are defined on an infinite 
strip (0, X) x (—oo, -boo) . We require that these nonperiodic functions be bounded on the strip 
0 < x < A, absolutely integrable and representable as a Fourier integral along y for any z, which, 
in particular, means that they vanish along any line x — const , u(x, y) — > 0 as y — > ±oo, 
u = (u, v , p, p). By virtue of boundary conditions (3.4), we expect that the solution of the dif- 
ference AP should converge on any rectangle (0, A) x (— y, y) to the corresponding fragment of 
such a solution to the nonhomogeneous counterpart of (3.2a) that satisfies (3.2b). (The RHS that 
would make (3.2a) nonhomogeneous is compactly supported.) Therefore, we can choose the values 
h x , h y , and Y such that the AP solution obtained in Dy for some compactly supported RHS will 
be arbitrarily close on (0, X) x (— y) D D xn to a function that can be (smoothly and uniquely) 
complemented on the entire R 2 to the vanishing at infinity solution of the inhomogeneous counter- 
part to (3.2a) that is generated by the same RHS. A detailed discussion of the type of convergence 
introduced above, as well as some estimates connecting h x , h v and T, are contained in [25]. 

In summary, we can say that instead of analytically transferring the boundary conditions from 
infinity to some closed finite boundary (e.g., the circle in (1.5) or the sphere in (2.3)), we now 
semi-analytically transfer boundary condition (3.2b) from infinity to dDy (see Figure 3.1) so that 
the solution to the AP on any fixed neighborhood of Di n can be as close to the solution of the 
inhomogeneous system (3.2a) (with the same compactly supported RHS) as desired. Clearly, con- 
sideration of the convergence only on some finite neighborhood of D{ n presents no loss of generality 
from the standpoint of constructing the DPM-based ABC’s. Indeed, this neighborhood can always 
be chosen to entirely contain r 1? and the BEP is always formulated with respect to a function 
concentrated near T. Therefore, the ABC’s obtained on the basis of the AP formulated above can 
be expected to be arbitrarily close to the exact ABC’s. In other words, these ABC’s meet the 
aforementioned fundamental requirement: one should be able to uniquely complement the solution 
calculated inside the finite computational domain to its infinite exterior so that the original problem 
is solved within the desired accuracy. 

The solution to the difference AP can be easily computed by means of the Fourier technique. 
The type of boundary conditions (3.4) (which are imposed separately for each wavenumber k ) makes 
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this choice of numerical method most appropriate. In so doing, one must solve system (3.3)-(3.4) 
for all k, k = We now briefly describe the fast algorithm for solving (3.3)— (3.4); this 

algorithm is delineated by Ryaben’kii and Tsynkov [91]. 

Denote by C the Unear space of eight- component complex- valued vectors; k £ C, 8 °m.k € C 
Vm, k. Conditions (3.4a) and (3.4b) select the subspaces C~ C C and C + C C, respectively; 
clearly, C~ H C + = 0, C = C~ © C + . Those solutions to the homogeneous counterpart of (3.3) 
that decrease as m — *• — oo belong to C~ at each node. Analogously, those solutions to the 
homogeneous counterpart of (3.3) that do not increase as m — * -f oo belong to C + at each node. 

Since the algorithm for solving (3.3)-(3.4) is the same for all k, k = —J, . . J, we drop the 
subscript k hereafter (as well as the “hat” ~ and the superscript 0, for simphcity). 

Let us now specify v l Q = 0 (the superscript l means left) and “integrate” (3.3) from left to right, 

(3-5a) v' l m = -Qv^ + A - 1 g m , m = l,...,M, 

implementing the projection 

(3.5b) vi, = vl - v'l, € C + , v'l € C~, 

onto the subspace C + at each step. Clearly, for any v'^ the representation v /l m = + v" l m (see 

(3.5b)) is unique. Analogously, let v r M = 0 (the superscript r means right) and “integrate” (3.3) 
from right to left, 


(3.6a) 


vl_ 1 = -Q- 1 v^ + B- 1 g m , 


m = M , . . . , 1 , 


(3.6b) 


m— 1 


— V 


IT 

m — 1 


— V 


tfT 

771 — 1 5 


i G 


C-, 


V 


nr 

771 — 1 


€ C + . 


Representation (3.6b) is also unique (for any v'^.j). 

The following proposition has been proven in [91]: 

Proposition 3.1. The vector-function v m d =l v r m + 0 < m < M, is the solution to (3.3) 

that satisfies (3.4 ) • 

The stability of the computational process (3.5)-(3.6) has been justified in [91] as well. It is 
also important to note that procedures (3.5) and (3.6) are both numerically cheap because they 
each require only two eighth-order matrix-vector multiplications at each step. The total cost is, 
therefore, O(M) operations for each k, k = — J, . . . , J. 

As soon as we formulate an AP with the desirable behavior of the solution in the far field 
and develop the numerical algorithm for solving this AP, further construction of the generalized 
potential, the BEP, and finally, the ABC’s, is straightforward. Basically, to obtain the discrete 
BEP the same steps that were described in Section 2 for the continuous Laplace equation should 
be repeated here for the specific finite-difference analogue to (3.2a). Then, the BEP is solved by 
means of some variational technique. In other words, for any distribution of u, v, p, and p along 
T, we find a continuation of these data to the grid boundary 7 that satisfies the BEP. Since the 
continuation from T to 7 is constructed on the basis of the Taylor expansion, then by solving the 
BEP we simply express normal derivatives of the solution on T in terms of boundary values of the 
solution itself. Finally, we use the grid density obtained (i.e., the solution to the BEP) to construct 
the generalized potential; the values of this potential are then interpolated from the Cartesian AP 
grid to T x . This process provides us with the desirable ABC’s, details of the procedure briefly 
described above can be found in [1, 25]. We only note that because the entire procedure is based 
on the linear formulation of the problem, the ABC’s can be finally obtained in a matrix form 
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(3.7) u = Tu , 

v 7 Iri lr’ 

which makes their implementation simple from an algorithmic standpoint and convenient for prac- 
tical computing. In (3.7), u represents the entire vector of u, v, p 7 and p along T, u represents 

r t I r i 

an analogous vector of unknowns along T i, and T is the operator (matrix) of boundary conditions 

that is calculated on the basis of the solution to the corresponding difference BEP. 

Let us now make a few remarks. First, formulation of the problem presented above deals with 
the true molecular Reynolds number and, therefore, does not account for turbulence. However, 
turbulent flows are undoubtedly most interesting for applications. Therefore, we qualitatively 
describe the far-field turbulent flow as a laminar flow of model fluid with a new “effective” viscosity 
(Boussinesq). An effective viscosity is chosen so that the corresponding model laminar solution 
coincides in the far field with an original turbulent wake-type solution. We also use Clauser’s 
conjecture to obtain a specific value of the constant in the expression for the effective viscosity and 
Tucker’s conjecture to account for the compressibility. A detailed description of this approximate 
approach for treating the turbulence in the far field can be found in the work of Tsynkov, Turkel, 
and Abarbanel [90]. Numerical results from [90] corroborate the statement that if applied in the far 
field the concept of effective turbulent viscosity can provide accurate computations without making 
the ABC’s more cumbersome. 

Second, the computation of boundary condition (3.7) (more precisely, the computation of 
matrix T from (3.7)) can be easily parallelized on a multi-processor machine. Indeed, the one- 
dimensional difference boundary-value problems (3,3)-(3.4) that correspond to different wavenum- 
bers k are obviously independent; therefore, they can be solved simultaneously on different pro- 
cessors of a parallel computer. This approach was practically implemented on an eight-processor 
CRAY Y-MP; such a parallelization appears to reduce the wall-clock time required for calculation 
of ABC (3.7) by approximately a factor of 5. 

Third, rather than using the uniform grid and the discrete Fourier transform with respect to j/, 
the stretched grid and, accordingly, the expansion with respect to a skew basis can be implemented. 
Numerical experiments (see below) indicate that with this approach high accuracy of the final results 
can still be maintained and, at the same time, the computer effort required to calculate the ABC’s 
can be drastically decreased. 

Referring the reader to the original papers [25], [1], [90], and [91] for more details regarding the 
construction of the DPM-based ABC’s for viscous flow computations, we now present some results of 
the numerical experiments with boundary conditions (3.7). We have combined the ABC (3.7) with a 
second-order finite- volume code [26, 27, 28] by Jameson, Schmidt, Turkel, and Swanson. The code is 
intended for calculating steady solutions to the Navier-Stokes equations; it uses multigrid iterations 
based on Runge-Kutta time stepping. The original (henceforth referred to as standard) treatment of 
the external boundary in the code [26, 27, 28] is local; this treatment is based on the extrapolation 
of characteristic and/or physical variables; the point-vortex correction can be employed as well. 
We have conducted a series of computations for different subsonic and transonic, laminar and 
turbulent flows around different airfoils. The results obtained in using the nonlocal ABC (3.7) were 
compared with those obtained on the basis of standard local boundary conditions with or without 
the point-vortex correction. The parameters of the specific computational variants, as well as the 
corresponding results, are reported in [1, 90]. The general conclusion that can be drawn from 
analyzing these results is that the nonlocal DPM-based ABC’s enable one to essentially shrink (by 
more than a factor of 10) the computational domain preserving the accuracy of computations or, 
conversely, to noticeably improve the accuracy for a fixed domain, i.e., keeping the computational 
cost at the same level. 
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In Table 3.1, we show some recent results that are not included in our previous work. Here, 
we compare the solutions (specifically, the dynamic lift C \ , dynamic drag Cd, and full drag Cd 
coefficients of a lifting airfoil under the non-zero angle of attack a) obtained for computational 
domains of different sizes and for different types of ABC’s, namely, nonlocal boundary conditions 
(3.7) and the point-vortex model. (For convenience, Table 3.1 presents both the actual values of 
Cu Cd, and Cd and the relative errors with respect to the corresponding asymptotic value; the 
latter is calculated for a domain with an “average radius” of 50 chords for each type of the ABC’s.) 

Table 3.1 

Comparison with point-vortex (p.-v.) model for RAE 2822 airfoil ; Mo = 0.73; Re o = 6.5 - 10 6 / a — 2.79°; basic 
grid 640 X 128 nodes; normal spacing 0.5 ■ 10“ 5 . 


Domain “radius” 

3 chords 

8 chords 

50 chords 

Grid 

600 x 104 

608 x 112 

640 x 128 

Type of ABC’ 

p.-v. 

(3.7) 

p.- V. 

(3.7) 

p.-v. 

(3.7) 

Cl 

0.8653 

0.8591 

0.8624 

0.8589 

0.8603 

0.8593 

relative error 

0.58% 

0.02% 

0.24% 

0.04% 

0% 

0% 

Cd x 10 

0.1203 

0.1263 

0.1209 

0.1261 

0.1255 

0.1260 

relative error 

4.14% 

0.24% 

3.67% 

0.08% 

0% 

0% 

Cd x 10 

0.1755 

0.1816 

0.1762 

0.1815 

0.1810 

0.1815 

relative error 

3.04% 

0.05% 

2.65% 

0% 

0% 

0% 


One can easily see that the asymptotic (50 chords) values of C\, Cd, and Cd that correspond to 
the different types of ABC’s are close to one another. However, as the artificial boundary approaches 
the airfoil, the discrepancy between the results produced by the two types of boundary conditions 
grows. We emphasize that the lift coefficients for both types of ABC’s are fairly stable with respect 
to the change in the size of the computational domain. As concerns the drag coefficient, it remains 
stable only for nonlocal ABC’s (3.7) (which are close to the exact ABC’s) and deteriorates for 
the point-vortex boundary conditions. This behavior is reasonable because lift actually drives the 
point- vortex model and drag is not taken into account by this model at all. 

Another important aspect of implementation of nonlocal ABC’s (3.7) is the influence they exert 
on the convergence rate of the multigrid iteration procedure. For some computational variants, the 
convergence for the DPM-based boundary conditions is up to 3 times faster than the convergence 
for a standard local technique. An example is shown in Figure 3.2, in which convergence dynam- 
ics relevant to the two types of ABC’s for the computation of a laminar subsonic flow past the 
NACA0012 airfoil are compared. Other examples can be found in [1, 90]. 

Finally, we note that boundary conditions (3.7) may also improve the robustness of the entire 
computational procedure. Indeed, the original iteration procedure [26, 27, 28] is rather sensitive 
to a number of factors that influence convergence (e.g., some geometric properties of the grid, as 
well as some iteration parameters). For certain computational variants, nonlocal ABC’s (3.7) can 
ensure good convergence when the standard procedure simply fails to converge. Again, specific 
examples can be found in [1, 90], 

We should also note that the computational cost of ABC’s (3.7) themselves is not high. Gener- 
ally, this expense never exceeds 10% of the overall cost of the computation. Clearly, in the case of 
strong convergence acceleration much can be gained from using the DPM-based ABC’s. However, 
even if the convergence acceleration is not as drastic as that shown in Figure 3.2, an additional 10% 
is still an acceptable cost increase if one takes into account the possibility to improve the accuracy 
and to use smaller computational domains. 

Summarizing our experience in developing and implementing nonlocal DPM-based ABC’s, we 
can say that the approach appears most promising in constructing effective and robust numerical 
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Fig. 3.2. Convergence history (log UpreitduailUoo versus number of multigrid cycles) for NACA0012 airfoil; 
Mo = 0.63; a = 2 0 ; Re = 5000; average radius of computational domain is approximately 6 chords. 


algorithms for solving external flow problems. The geometric universality of this method along 
with the algorithmic simplicity of its implementation on one hand, and the high accuracy of the 
obtained results on the other hand can make it an effective tool in modern CFD. 

There are several possible extensions of the described approach. The one most demanded 
by current computational practice is an extension to 3D steady-state flows. Although relevant 
experiments have not yet been conducted, such an extension seems feasible. The general framework 
of the algorithm would remain the same. As concerns the technical portion, it would certainly 
become more cumbersome in particular because of the need to consider 2D surfaces rather than 
ID curves as artificial boundaries. However, highly accurate and robust ABC’s may be even more 
crucial for computations in three dimensions than in two dimensions since in three dimensions there 
are, generally speaking, no simple models, such as the point-vortex model, that can noticeably 
improve the results provided by standard local techniques. 

Another interesting area is the combined implementation of nonlocal ABC’s and multigrid 
iteration procedures. As mentioned above, our computations [1] and the computations by Ferm 
[24] confirm that the nonlocal ABC’s can be particularly effective if used in combination with 
multigrid. However, many modern multigrid techniques in CFD are often not optimal themselves. 
A massive effort is currently underway toward the construction of new discretizations that may 
improve the efficiency of multigrid solvers. Development of highly accurate boundary conditions 
may essentially contribute to achieving this goal. 

However, the most challenging problem in our opinion is the extension of the DPM-based 
approach to the treatment of unsteady flows. A particular class of such problems, specifically, 
the flows that oscillate in time, has been studied by Tsynkov [92]. In practice, this formulation 
originates from the well-known problem of an oscillating airfoil or from the problem of the time- 
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periodic injection of fluid into the boundary layer, which according to the experimental data of 
Seifert et al. [93] may essentially improve the performance of airfoil. In [92], we linearize the 
governing equations in the far field and then, since the main frequency is known, implement the 
Fourier transform in time. If the problem is already discretized in time, then we end up with a 
finite family of steady-state problems (generally speaking, with complex coefficients); each of these 
problems can be treated in a way analogous to the one described above. Additionally, we provide 
in [92] a more thorough study of the solvability of AP and propose a new technique for practically 
setting the DPM-based ABC’s; this new technique does not require solution of the BEP. Obviously, 
the nonlocal character of the ABC’s [92] in time is restricted by the value of one period. Moreover, 
the theoretical analysis in [92] is in a sense less cumbersome than the analysis in [17] because at 
the end we deal with stationary rather than general time-dependent equations. Therefore, we can 
say that the time-periodic case occupies an intermediate position between steady-state and truly 
time-dependent problems. In regard to the latter, the primary difficulty is how to effectively “cut 
off the tail” in time. This issue is the subject of a future investigation. 
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